
function H = HessLL_cfs(theta,~,X_PVEM,X_PRI,ES0,ES1,ES2)

P1 = exp(X_PVEM * theta + ES1 - ES0);
P2 = exp(X_PRI * theta + ES2 - ES0);
P0 = 1 + P1 + P2;
P1 = P1 ./ P0;
P2 = P2 ./ P0;

H = - X_PVEM' * diag(P1.*(1-P1)) * X_PVEM + ...
    X_PVEM' * diag(P1.*P2) * X_PRI + ...
    X_PRI' * diag(P1.*P2) * X_PVEM + ...
    - X_PRI' * diag(P2.*(1-P2)) * X_PRI;

H=-H;



